Species Abundance Patterns in Complex Evolutionary Dynamics 
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An analytic theory of species abundance patterns (SAPs) in biological networks is presented. 
The theory is based on multispecies replicator dynamics equivalent to the Lotka-Volterra equation, 
with diverse interspecies interactions. Various SAPs observed in nature are derived from a single 
parameter. The abundance distribution is formed like a widely observed left-skewed lognormal 
distribution. As the model has a general form, the result can be applied to similar patterns in other 
complex biological networks, e.g. gene expression. 

PACS numbers: 87.23.-n,75.10.Nr,87.10.+e,87.90.+y 



If we investigate the number and populations of species 
in an ecosystem, we can observe universal character- 
istic patterns in that ecosystem. How to clarify the 
mechanisms underlying those species abundance patterns 
(SAPs) has been one of the 'unanswered questions in ecol- 
ogy in the last century [l|' even though the knowledge 
obtained from it would affect vast areas of nature con- 
servation. Various models have been applied to ecosys- 
tem communities where species compete for niches on a 
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Hq , but these models have left the more complex systems 
a mystery. Such systems occur on multiple trophic lev- 
els and include various types of interspecies interactions, 
such as prey-predator relationships, mutualism, compe- 
tition, and detritus food chains. Although SAPs are ob- 
served universally in nature, their essential parameters 
have not been fully clarified. 

I consider, then, a widely adopted model of biological 
networks represented by the so-called V-species replica- 
tor equation (RE) |l7j : 
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to calculate the abundance Xi(t)(& [0, N}) of species 
i(= 1,2 ... ,N). Here we assume that (Jij) is a time- 
independent random symmetric ( Jij = Jji) matrix whose 
elements have a normal distribution with mean m(> 0) 
and variance J 2 /N as 
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Self-interactions are all set to a negative constant as 
Ju = —u(< 0). Note that the essential parameter is 
unique as p = (u + m)/J because the transformation of 
the interaction K^j = (Jy — m)/J does not change the 
trajectory of the dynamics (J). Although ecologists do 
not generally believe in the randomness of interspecies 
interactions in nature, the discipline has been affected 
by the random interaction model [l8| as a prototype of 
complex systems. 



The RE appears in various fields [11]]. In sociobiol- 
ogy, it is a game dynamical equation for the evolution of 
behavioral phenotypes; in macromolccular evolution, it 
is the basis of autocatalytic reaction networks (hypercy- 
cles) ; and in population genetics it is the continuous-time 
selection equation in the symmetric (Jy = Jji) case. The 
symmetric RE also corresponds to a classical model of 
competitive community for resources E^- 

Particularly in the context of ecology, the N—l species 
Lotka-Volterra (LV) equation 
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is equivalent to the N species RE 17] . That is, the abun- 
dance yi and the parameters in the corresponding LV are 
described by those in the present RE model as, 
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where the 'resource' species M (yM = 1) can be arbi- 
trarily chosen from N species in the RE. The ecological 
interspecies interactions {bij)(i ^ j) have a normal dis- 
tribution with mean and variance 2J 2 /N from Eq. 10, 
and they are no longer symmetric (bij ^ bji). The 
present model therefore describes an ecological commu- 
nity with complex prey-predator interactions ((bij , bji) — > 
(+, — ) or (— , +)), mutualism (+,+) and competition 
(—,—). Moreover, a community can have a 'loop' 
(detritus) food chain ((bij, bji) — * (+,-), (b jk ,b kj ) — > 
(+,—), (pki-, bik) — > (+, ^ ))• The intraspecific interac- 
tion bu turns out to be related to the intrinsic growth 
rate r% as bu = Ju — Jmi — —u — JiM = ~ r i and is 
therefore competitive (bu < 0) for producers (r, > 0) or 
mutualistic (bu > 0) for consumers (r^ < 0). 

By Eq. (J5J, the intrinsic growth rates also have a nor- 
mal distribution with mean u + m and variance J 2 /N. 
The probability at which r.; is positive-that is, that the 



2 



i-th species is a producer-is therefore given by the error 
function, 

r°° dt 

Prob(r 4 > 0) = / -= cxp (-t 2 ) (7) 
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Consequently, the parameter p can be termed as the 'pro- 
ductivity' of a community because the larger the p, the 
greater the number of producers. The parameter p is 
also connected to the maturity of an ecosystem because 
m increases in time in an evolutionary model |20| . 

The symmetry (J.^ = Jji) makes the average fitness 

/ = Ylj k JjkXjXk (the second term of the r.h.s. of 
Eq. a Lyapunov function |17|. which is a nondc- 
creasing function of time in dynamics Therefore, 
every initial state converges to a local maximum of / as 
t — > oo. Interpreting 3~C = —\j as an energy function, we 
can study macroscopic functions like free energy at such a 
maximum by using the tec hniq ue of statistical mechanics 
of random systems [U |H Hi lU H| . 

Information on equilibrium states of dynamics Q at 
t — > oo is derived from the zero-temperature limit of free 
energy density, 

/ = -Jim Jim -^- 3 [lnZ] J (8) 
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where [• -Jr = J^° OQ dJijP(Jij)(. . .) denotes the 'sample 
average |2l|' over random interactions. The Dirac delta 
function in Eq. @ reflects the conservation of total abun- 
dance Xi(t) — N satisfied at any t in Eq. The 
calculation of Eq. JSJ is similar to calculations in the 
previous works jU 0, H^, |2(| and yields mean field 
equations for the order parameters q and v as 
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where A = ^/q(p — 2v) and Dz = dzexp(— z 2 /2)/v / 2~7r. 
The resulting equations turn out to be formally the same 
as the case where m — and J = 1 j2^|. For each value 
of p, Eqs. 1)1 and (|ll)l are solved numerically. 

Among macroscopic functions calculated in the present 
framework, the most significant for a theory of SAPs 
is the survival function a p (x) = jj^2f &{xi — x), the 
proportion of species whose abundance is larger than x, 
where 6(z){= l(z > 0);0(z < 0)) is the step function. 
Similar to the fashion in which free energy was calculated, 
the survival function a p {x) is analytically calculated and 
represented by the order parameters as 
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FIG. 1: Diversity a v (x — 0, 1) as a function of p of log-log 
scales. Black circles show numerical solutions of a p (l) aver- 
aged over 50 samples of (Jij) for Eq. Q with N = 2048 and 
p = 0.1, 0.2, 0.3, 0.4, 0.5, V2/2, 1, \/2, 2, 3. Error bars indicate 
the maximum and minimum values found in the samples. 
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where the definition of the step function above is given 
by 6{z) = l(z > 0);0(z < 0)). 

The resulting function a p (0) = v(p — v) and a p (l) of p 
can be termed 'diversity', i.e., the proportion of nonex- 
tinct species and that of the species with abundance 
larger than unity, respectively, as depicted in Fig. ^ 
This demonstrates a typical positive correlation between 
productivity and diversity [27| . Numerical results for 
a p (l) are also depicted in Fig. ^ f° r comparison. We 
see good agreement between the analytical and the nu- 
merical results for p > 1, while some deviations appear 
for small values of p. This small- value deviation is at- 
tributable to the occurrence of replica symmetry break- 
ing (RSB)|11|23 for p < V2, which yields a number of 
metastable states of Eq. (JHJ), and the replicator dynamics 
Q essentially converges to not only a ground state of 
JHJ) but also to the metastable states. Since the energy 
IK and the diversity are both nonincreasing functions of 
time in dynamics the mean- field results here give a 
lower minimum of diversity. Interestingly, the metastable 
states enhance the diversity. The analysis of RSB is ex- 
pected to improve the quantitative agreement |23j . 

Note that a p (x) is also represented as a function of 
species rank n: 



a p (x) 
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(n = 1,2, . . . , S < N) if the species abundance is ranked 
in descending order, as in a;'- 1 -' > x^> > ■ ■ ■ > x^ > 
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FIG. 2: Rank-abundance relations as a function of produc- 
tivity p on normal-log scales. 



• • • > > 0. As the function a p (x) is a nonincreasing 
monotonic function, the species abundance relation, i.e., 
the abundance i'"' as a function of a rank n, is given 
by the inverse function of a p (x) as x^ n ' = x p (n/N) = 
a p (x), depicted in Fig.|2lfor some values of p. We ob- 
serve two typical SAPs in different regions an d with 
different species compositions |6J: one is a straight line 
like the geometric series for a small value of p, and 
the other consists of sigmoid curves on a logarithmic ver- 
tical axis for some range of p. This latter SAP denotes a 
lognormal-like abundance distribution. Remarkably, the 
transition of the SAPs from low p to high is identical 
to the observed transition from low- to high-productivity 
areas; that is, from a species-poor area such as an alpine 
or polar region to a species-rich tropical rain forest [12j . 
The transition also corresponds to the secular variation 
of SAPs observed in abandoned cultivated land Q ■ This 
supports the contention that p (orm) is a maturity pa- 
rameter, as is suggested by an evolutionary model 12(1 ]. 

The abundance distribution is also derived from the 
survival function. As C p (x) = 1 — a p (x) is a cumulative 
distribution function of abundance, the abundance dis- 
tribution is given by the derivative F p (x) = dC p (x)/dx 
and 
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where the second term denotes the rate of species ex- 
tinction. The first term is a normal distribution but 
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FIG. 3: Abundance distribution per 'natural' octave ln(a;). 
Functions F p (x)x(d(\n(x)) = F p (x)dx) for some values of p 
are depicted, whereas Preston originally defined octaves as 
logarithms to base 2 



not a lognormal distribution. Nevertheless, the curves 
in Fig. |2l demonstrate a typical sigmoid pattern on a log- 
arithmic vertical axis. This pattern indicates the coex- 
istence of very abundant species with rare ones. This 
multiscale of abundance is intuitively understood by a 
divergent behavior of the variance a 2 = q/(p — v) 2 of 
F p (x) for small p because q — > oo and v — > for p — > 0. 
Moreover, the mode of F p (x) per 'natural' octave 
ln(x) is always a positive value (as shown in Fig. I^J) at 
x* = f (A + \/A 2 -I- 4) > 0, which denotes a unimodal 
distribution. Indeed, the mode diverges as 
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for p — > 0. As a result, the abundance distribution is a 
truncated normal distribution with a large variance a 2 — > 
oo and a negatively divergent mean /i = 9 ^Z„"^ ^ — oo 
satisfying - = — > for p — > 0. This is why the abun- 
dance distribution per octave looks like a left-skewed log- 
normal distribution [13 in Fig. El 

Moreover, we derive an analytical expression for 
the population x max of the most abundant species, 
the position of the individual curve mode xn, and 
therefore the ratio of their logarithm defined as 7 = 
log(xN)J\og(x max ). According to the canonical hypoth- 
esis |5j, |8( , the parameter 7 takes a value near unity in 
various real communities. To check the validity of the 
canonical hypothesis, we first need to evaluate an ex- 
pected value of the most abundant species x max . From 
the definition of x max , that is N F p (x max ) — 1, and the 
conservation of the total abundance J F p (x)xdx = N, 
which is equivalent to Xi — N, we obtain 



q(p -v) + f7 y 2 ln(^^ »- +Aa '- (0) ) 
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On the other hand, the mode of the individual curve 
F p (x)x per octave is given by xm = § (A + %/A 2 + 8), 
and finally, the parameter 7 = log(a;jv)/log(a; maa ;) is 
evaluated by substituting the values of the order pa- 
rameters q and v for each value of p. In the present 
model, 7 is a monotonically increasing function of p and 
0.96 < 7 < 1.04 for 0.1 < p < 0.6, denoting that the 
canonical hypothesis is supported in the range of p giv- 
ing the typical SAPs in Fig. |3 Although the canonical 
hypothesis was demonstrated to be merely a mathemati- 
cal consequence of lognormal distribution 8] rather than 
anything biological, it is noteworthy that the lognormal- 
like abundance distribution with 7 ~ 1 derives from basic 
ecological dynamics. This still suggests a biological foun- 
dation for the hypothesis in a large complex ecosystem, 
in the same way that a biological foundation was indi- 
cated for the theory of a local competitive community 
0. 

In the present model, all species coexist only in the 
limit p — > 00, that is, in the trivial cases in which in- 
terspecies interactions are negligible (J « u ) or ho- 
mogeneous (J — > 0), thereby giving aoo(x) = 9(1 — x), 
= XQO (n/N) = 1 for all n and i^x) = S(x - 1). 

The present theory seeks to capture the influence of 
productivity on the SAPs under the assumption that 
all species interact randomly; nevertheless, this assump- 
tion itself is never justified because it ignores a biological 
correlation between interactions produced by evolution. 
However, note that the randomness is assumed only for 
an initial state with N species in Eq. (Q. Actually, the 
simulation reveals the resulting interactions of nonextinct 
species to be nonrandom: every sample for p € [0.1, 3] in 
Fig. n evolves to only flora, Vi r, > 0. Moreover, by 
ordering species as r< > Tj for any i < j, we observe a 
hierarchy: there are only three types of interactions, that 
is, mutualism (bij,bji) = (+,+), competition (— ,— ) and 
exploitation of i on j as (+,—), but no reverse (—,+). 
This suggests the applicability of the present model to a 
plant community. 

It has been demonstrated that empirically supported 
patterns are derived from a single parameter of general 
population dynamics. This not only suggests the im- 
portance of globally coupled biological interactions in a 
large assemblage but also provides a unified viewpoint 
on mechanisms of similar patterns observed in other bio- 
logical networks with complex interactions; for example, 
a lo gno rmal abundance distribution of a protein in cells 
[jjl I29L |30| . which is revealed by gene expression net- 
works. 
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